Set up

Install and load required packages

devtools::install_github("CWWhitney/uncertainty")
library(DiagrammeR)
library(decisionSupport)
library(ggplot2)
library(kableExtra)
library(magrittr)
library(plyr)
library(rmarkdown)
library(uncertainty)
library(chillR)
library(raster)

Load emission calculation function

The formulas come from Chapter 10 of the 2006 IPCC Guidelines for National Greenhouse Gas Inventories (Dong et al. 2006) (Volume 4 - Agriculture, Forestry and Other Land Use). The document can be found on the IPCC website. The chapter is part of the book Good Practice Guidance and Uncertainty Management in National Greenhouse Gas Inventories (IPCC 2021), also on the the IPCC website. Many variables are based on the work done in Variation in the carbon footprint of milk production on smallholder dairy farms in central Kenya (Wilkes et al. 2020) and in Methods and guidance to support MRV of livestock emissions (Wilkes et al. 2019). We run the model with the decisionSupport package (Luedeling et al. 2022).

# Generate the GHG emissions function 
source("GHG_emissions_function.R", local = knitr::knit_global())

This equation takes a wide range of inputs. Many of these are usually not known precisely, yet the classic IPCC equation does not account for potential errors and uncertainties. Here, we express this uncertainty by defining distributions for each variable. The following table shows all the variables, as well as examples of distributions that may describe their quantities. Distributions are defined according to the distribution shape, as well as upper and lower bounds of 90% confidence intervals.

### # This project will be stored in relative file paths (i.e. same or lower folders)

 input_table <- "data/livestock_ghg_input_table.csv"

options(knitr.kable.NA = '', encoding = "Windows-1252")
 read.csv(input_table)[,c(1:2,4:6)]   %>%
    kableExtra::kbl(digits = 2)  %>%
    kableExtra::kable_classic_2(full_width = F)
Description lower upper distribution variable
Feeding system and corresponding activity coefficients
Feeding system (1 Stall - 2 Pasture - 3 Grazing large area) 1.00 1.00 const feeding_system
Acticity coefficient Ca - Stall feeding 0.00 0.00 const Ca_stall
Acticity coefficient Ca -Pasture feeding 0.14 0.20 posnorm Ca_pasture
Acticity coefficient Ca - Grazing large area 0.31 0.41 posnorm Ca_large_grazing
Herd composition
Number of ‘intact’ mature bulls (>36 months) (animal type 1) 5.00 5.00 const N_animal_type_1
Number of mature castrates (>36 months) (animal type 2) 2.00 2.00 const N_animal_type_2
Number of males 12-36 months (animal type 3) 4.00 4.00 const N_animal_type_3
Number of lactating cows (animal type 4) 5.00 5.00 const N_animal_type_4
Number of lactating cows (animal type 5) 3.00 3.00 const N_animal_type_5
Number of lactating cows (animal type 6) 2.00 2.00 const N_animal_type_6
Number of weaned young females <=12 (animal type 7) 4.00 4.00 const N_animal_type_7
Number of weaned young males <=12 (animal type 8) 3.00 3.00 const N_animal_type_8
Number of females 12-36 months (animal type 9) 5.00 5.00 const N_animal_type_9
Number of pre-weaned female calves (animal type 10) 6.00 6.00 const N_animal_type_10
Number of pre-weaned male calves (animal type 10) 4.00 4.00 const N_animal_type_11
Number of first-time pregnant cows (animal type 12) 3.00 3.00 const N_animal_type_12
Sex-specific and (maintenance) coefficients
Sex-specific coefficient - intact males 1.02 1.38 posnorm C_intact_males
Sex-specific coefficient - castrates 0.85 1.15 posnorm C_castrates
Sex-specific coefficient - females 0.68 0.92 posnorm C_females
Maintenance coefficient (Cfi) - Cattle/Buffalo (non-lactating cows, steers and juveniles) 0.27 0.37 posnorm maintenance_coefficient_Cfi_other
Maintenance coefficient (Cfi) - Cattle/Buffalo (lactating cows) 0.33 0.44 posnorm maintenance_coefficient_Cfi_lact_cow
Maintenance coefficient (Cfi) - Cattle/Buffalo (bulls) 0.31 0.43 posnorm maintenance_coefficient_Cfi_bulls
Mature body weight females 279.20 418.80 posnorm mature_weight_MW_female
Mature body weight males 386.40 579.60 posnorm mature_weight_MW_male
Weight (gain) data by animal type
Mean live weight of animal type 1 99.50 537.68 posnorm live_weight_LW_1
Mean live weight of animal type 2 84.40 468.84 posnorm live_weight_LW_2
Mean live weight of animal type 3 81.50 374.07 posnorm live_weight_LW_3
Mean live weight of animal type 4 69.20 480.83 posnorm live_weight_LW_4
Mean live weight of animal type 5 59.20 454.38 posnorm live_weight_LW_5
Mean live weight of animal type 6 62.40 488.65 posnorm live_weight_LW_6
Mean live weight of animal type 7 61.90 233.83 posnorm live_weight_LW_7
Mean live weight of animal type 8 45.20 191.35 posnorm live_weight_LW_8
Mean live weight of animal type 9 90.80 402.37 posnorm live_weight_LW_9
Mean live weight of animal type 10 23.60 93.12 posnorm live_weight_LW_10
Mean live weight of animal type 11 18.90 96.89 posnorm live_weight_LW_11
Mean live weight of animal type 12 63.90 425.82 posnorm live_weight_LW_12
Average daily weight gain (kg/day) of animal type 1 0.00 0.00 const weight_gain_WG_1
Average daily weight gain (kg/day) of animal type 2 0.00 0.00 const weight_gain_WG_2
Average daily weight gain (kg/day) of animal type 3 0.12 0.18 posnorm weight_gain_WG_3
Average daily weight gain (kg/day) of animal type 4 0.00 0.00 const weight_gain_WG_4
Average daily weight gain (kg/day) of animal type 5 0.00 0.00 const weight_gain_WG_5
Average daily weight gain (kg/day) of animal type 6 0.00 0.00 const weight_gain_WG_6
Average daily weight gain (kg/day) of animal type 7 0.34 0.50 posnorm weight_gain_WG_7
Average daily weight gain (kg/day) of animal type 8 0.28 0.40 posnorm weight_gain_WG_8
Average daily weight gain (kg/day) of animal type 9 0.21 0.31 posnorm weight_gain_WG_9
Average daily weight gain (kg/day) of animal type 10 0.34 0.50 posnorm weight_gain_WG_10
Average daily weight gain (kg/day) of animal type 11 0.28 0.40 posnorm weight_gain_WG_11
Average daily weight gain (kg/day) of animal type 12 0.21 0.31 posnorm weight_gain_WG_12
Milk yield; fat content; pregnancy coefficients
Average milk yield (kg/day) of animal type 4 5.11 8.35 posnorm milk_yield_4
Average milk yield (kg/day) of animal type 5 5.11 8.35 posnorm milk_yield_5
Average milk yield (kg/day) of animal type 6 5.11 8.35 posnorm milk_yield_6
Milk fat (%) 3.60 4.40 posnorm milk_fat
Maximum pregnancy coefficient Cp (scaled to year) 0.10 0.10 const Cp_0
Pregnancy coefficient for animal type 4 (share of full coefficient) 0.00 0.89 posnorm Cp_rate_4
Pregnancy coefficient for animal type 5 (share of full coefficient) 0.11 1.00 posnorm Cp_rate_5
Pregnancy coefficient for animal type 6 (share of full coefficient) 0.01 0.99 posnorm Cp_rate_6
Emissions from externally produced feed
Feed production 834.56 834.56 const feed_prod_CO2
Feed transport 5.20 5.20 const feed_trans_CO2
Conversion parameters
Digestible energy expressed as a percentage of gross energyy (De) 58.67 59.04 posnorm DE_perc
Percent crude protein in diet 5.00 15.00 posnorm perc_crude_protein_diet
Methane conversion rate (0-1) 0.05 0.08 posnorm Y_m
emission factor for N2O emissions from nitrogen leaching and runoff, kg N2O-N/kg N leached and runoff 0.01 0.01 const EF5
emission factor for N2O emissions from atmospheric deposition of nitrogen on soils and water surfaces 0.01 0.01 const EF4
urinary energy expressed as fraction of Gross Energy intake 0.04 0.04 const urinary_energy
the ash content of manure calculated as a fraction of the dry matter feed intake 0.08 0.08 const ash_content
Constants for methane production and conversion
Maximum methane producing capacity of the manure (animal type 1) 0.11 0.15 posnorm Bo_animaltype_1
Maximum methane producing capacity of the manure (animal type 2) 0.11 0.15 posnorm Bo_animaltype_2
Maximum methane producing capacity of the manure (animal type 3) 0.11 0.15 posnorm Bo_animaltype_3
Maximum methane producing capacity of the manure (animal type 4) 0.11 0.15 posnorm Bo_animaltype_4
Maximum methane producing capacity of the manure (animal type 5) 0.11 0.15 posnorm Bo_animaltype_5
Maximum methane producing capacity of the manure (animal type 6) 0.11 0.15 posnorm Bo_animaltype_6
Maximum methane producing capacity of the manure (animal type 7) 0.11 0.15 posnorm Bo_animaltype_7
Maximum methane producing capacity of the manure (animal type 8) 0.11 0.15 posnorm Bo_animaltype_8
Maximum methane producing capacity of the manure (animal type 9) 0.11 0.15 posnorm Bo_animaltype_9
Maximum methane producing capacity of the manure (animal type 10) 0.11 0.15 posnorm Bo_animaltype_10
Maximum methane producing capacity of the manure (animal type 11) 0.11 0.15 posnorm Bo_animaltype_11
Maximum methane producing capacity of the manure (animal type 12) 0.11 0.15 posnorm Bo_animaltype_12
Methane conversion factor - pasture 0.00 0.01 posnorm MCFpasture
Methane conversion factor - spreading 0.00 0.01 posnorm MCFspread
Methane conversion factor - drylot 0.01 0.03 posnorm MCFdrylot
Methane conversion factor - solid storage 0.03 0.07 posnorm MCFsolidstore
Methane conversion factor - composted 0.01 0.04 posnorm MCFcomposted
Methane conversion factor - liquid 0.40 1.20 posnorm MCFliquid
Methane conversion factor - biogas 0.05 0.16 posnorm MCFbiogas
Methane conversion factor - burning 0.05 0.15 posnorm MCFburn
Methane conversion factor - sold 0.00 0.00 const MCFsold
Decision tree for manure management (no longer needed)
question 1: how much of the manure is centrally collected? 0.10 0.90 tnorm_0_1 share_centrally
question 2: of manure that’s not centrally collected, how much is burned? 0.10 0.90 tnorm_0_1 share_central_burn
question 3: of centrally collected manure, how much is spread daily on cropland or pasture? 0.10 0.90 tnorm_0_1 share_spread
question 4: how much of the manure that’s not spread is sold? 0.10 0.90 tnorm_0_1 share_sold
question 5: how much of the stored manure that’s not sold is stored in liquid form? 0.10 0.90 tnorm_0_1 share_liquid
question 6: how much of the liquid manure is converted to biogas? 0.10 0.90 tnorm_0_1 share_biogas
question 7: how much of the solid manure is composted? 0.10 0.90 tnorm_0_1 share_compost
question 8: how much of the solid manure is stored in dry lots? 0.10 0.90 tnorm_0_1 share_drylot
Manure management system
Share of manure management (pasture) 0.30 0.30 const mms_pasture
Share of manure management (spread) 0.10 0.10 const mms_spread
Share of manure management (dry lot) 0.10 0.10 const mms_drylot
Share of manure management (solid storage) 0.10 0.10 const mms_solidstore
Share of manure management (composted) 0.10 0.10 const mms_composted
Share of manure management (liquid) 0.10 0.10 const mms_liquid
Share of manure management (biogas) 0.10 0.10 const mms_biogas
Share of manure management (burned) 0.10 0.10 const mms_burn
Share of manure management (burned) 0.00 0.00 const mms_sold
Manure emission coefficients by management
Manure nitrogen direct emission coefficients (Table 10.21 from the IPCC guidelines)
Manure nitrogen direct emission coefficient (pasture) 0.00 0.04 posnorm mndec_pasture
Manure nitrogen direct emission coefficient (spread) 0.00 0.00 const mndec_spread
Manure nitrogen direct emission coefficient (dry lot) 0.00 0.04 posnorm mndec_drylot
Manure nitrogen direct emission coefficient (solid storage) 0.00 0.02 posnorm mndec_solidstore
Manure nitrogen direct emission coefficient (composted) 0.00 0.02 posnorm mndec_composted
Manure nitrogen direct emission coefficient (liquid) 0.00 0.02 posnorm mndec_liquid
Manure nitrogen direct emission coefficient (biogas) 0.00 0.00 const mndec_biogas
Manure nitrogen direct emission coefficient (burned) 0.00 0.00 const mndec_burn
Manure nitrogen direct emission coefficient (sold) 0.00 0.00 const mndec_sold
Manure nitrogen volatilization emission coefficients (Table 10.22 from the IPCC guidelines)
Manure nitrogen volatilization emission coefficient (pasture) 0.02 0.02 const mnvc_pasture
Manure nitrogen volatilization emission coefficient (spread) 0.05 0.60 posnorm mnvc_spread
Manure nitrogen volatilization emission coefficient (dry lot) 0.20 0.50 posnorm mnvc_drylot
Manure nitrogen volatilization emission coefficient (solid storage) 0.10 0.65 posnorm mnvc_solidstore
Manure nitrogen volatilization emission coefficient (composted) 0.14 0.70 posnorm mnvc_composted
Manure nitrogen volatilization emission coefficient (liquid) 0.09 0.36 posnorm mnvc_liquid
Manure nitrogen volatilization emission coefficient (biogas) 0.00 0.00 const mnvc_biogas
Manure nitrogen volatilization emission coefficient (burned) 0.00 0.00 const mnvc_burn
Manure nitrogen volatilization emission coefficient (sold) 0.00 0.00 const mnvc_sold
Fraction of managed manure nitrogen losses for livestock category T due to runoff and leaching during solid and liquid storage of manure (typical range 0.01-0.20)
Fraction of managed manure nitrogen losses … (pasture) 0.01 0.73 posnorm frac_leach_pasture
Fraction of managed manure nitrogen losses … (spread) 0.01 0.73 posnorm frac_leach_spread
Fraction of managed manure nitrogen losses … (dry lot) 0.01 0.73 posnorm frac_leach_drylot
Fraction of managed manure nitrogen losses … (solid storage) 0.01 0.73 posnorm frac_leach_solidstore
Fraction of managed manure nitrogen losses … (composted) 0.01 0.73 posnorm frac_leach_composted
Fraction of managed manure nitrogen losses … (liquid) 0.01 0.73 posnorm frac_leach_liquid
Fraction of managed manure nitrogen losses … (biogas) 0.00 0.00 const frac_leach_biogas
Fraction of managed manure nitrogen losses … (burned) 0.00 0.00 const frac_leach_burn
Fraction of managed manure nitrogen losses … (sold) 0.00 0.00 const frac_leach_sold

Monte Carlo simulation

With the model equation and the input table defined above, we can run a Monte Carlo simulation that returns a probability distribution for the greenhouse gas emissions resulting from the particular herd that is specified in the input table. We use the functions provided by the decisionSupport package (Luedeling et al. 2022) in the R programming language (R Core Team 2022).

GHG_emissions<-mcSimulation(estimate=estimate_read_csv(input_table),
                            model_function=ghg_emissions,
                            numberOfModelRuns=1000,
                            functionSyntax="plainNames")

Running a Monte Carlo simulation with scenarios

The mcSimulation function takes random draws for each of the input variables from the distributions specified in the input table. The distributions that are used are the same for every run of the Monte Carlo simulation. This may be desirable in homogeneous populations where all combinations of random values are plausible, but it often misses the mark where populations are composed of particular member types with specific characteristics. To account for populations that are heterogeneous, which is the case for livestock-keeping households in Kenya, we use the scenario_mc function, which allows specifying particular subgroups. Defining such scenarios can be desirable when it’s necessary to simulate outcomes for particular farm types, climate scenarios etc. To do this, we would normally have to define separate input tables and run separate simulations. This is quite inconvenient. Here, we use a function that can do this job automatically.

The scenario_mc function in the decisionSupport package (Luedeling et al. 2022) consists of a wrapper around the mcSimulation function. The function is able to modify the input table according to a scenario file we provide. This file contains information on which variables should be modified for each scenario and how many model runs should be included in the simulation. The function is essentially the same as the mcSimulation function, but with an additional scenarios option which imports the contents of a .csv file.

The .csv file for the scenarios option specifies variables that should be modified for each of the scenarios, as well as the number of model runs for it. Each column beyond the first two columns in the file contain information for a particular scenario. This is restricted to the parameters that should be changed for this scenario.

Here’s an example scenarios file:

knitr::kable(read.table("data/scenarios.csv", sep = ",", fileEncoding = "UTF-8-BOM",header=TRUE))
Variable param Scen1 Scen2
Runs x 100 300
N_animal_type_1 both 100 4000
N_animal_type_2 lower 100 10
N_animal_type_2 upper 1000 20
N_animal_type_2 distribution posnorm norm

The first row (“Runs” in the first column) indicates the number of runs. Each subsequent row indicates which value from the input table should be modified, so the string in the first column has to correspond to a variable name in the input table. The second column indicates which values should be changed, with options being “lower,” “upper,” “distribution” and “both.” If “both” is selected, both the lower and upper bounds are modified (this should only be done for variables witha constant distribution, since the rplacement will be for both the upper and lower bound).

We define subgroups based on the populations of households described by a household survey in Kenya (Wilkes et al. 2020, 2019). Each of the farms encountered in this survey were converted into one scenario, for which Monte Carlo simulation is then used to produce multiple replicates.

Household-based scenarios

The information for the household is extracted from a spreadsheet that contains all the survey responses (Wilkes et al. 2020, 2019). Extracting all relevant information from these spreadsheets is a slightly complex task, which is accomplished by the following code.

# import the data from the survey and ensure that the .csv encoding can be read by R

Kenyaherds<-read.table("data/Kenya_baseline_individual_cow_wNotes.csv", sep = ",",
                       fileEncoding = "UTF-8-BOM", header=TRUE)

# count the number of animals for each type

animals<-table(Kenyaherds[,c("QQNO","animal_type")])

variables_to_update_const<-c(paste0("N_animal_type_",1:12),
                             "feeding_system")

scenarios<-data.frame(Variable=variables_to_update_const,param="both")

variables_to_update_posnorm<-c(paste0("milk_yield_",4:6),
                             paste0("Cp_rate_",4:6),
                             "feed_prod_CO2",
                             "feed_trans_CO2",
                             paste0("live_weight_LW_",1:12),
                             paste0("weight_gain_WG_",1:12),
                             "DE_perc",
                             "perc_crude_protein_diet")

scenarios<-rbind(scenarios,data.frame(Variable=rep(variables_to_update_posnorm,each=3),param=c("lower","upper","distribution")))


for(i in 1:nrow(animals))
  {coln<-paste0("Farm_",row.names(animals)[i])
   scenarios[,coln]<-NA
     scenarios[which(scenarios$Variable %in% paste0("N_animal_type_",1:12)),coln]<-animals[i,]
}

# for some animal-scale variables (e.g. milk yield), we need averages per farm

### error estimates for variables
milk_yield_error<-0.275 # Migose et al. (2020, https://www.sciencedirect.com/science/article/abs/pii/S1871141319307371 )
feed_trans_error<-0.25 # Vellinga et al (2013, https://edepot.wur.nl/254098)
feed_prod_error<-0.25 # Vellinga et al (2013, https://edepot.wur.nl/254098)
live_weight_error<-0.145 # Goopy et al. (2018, https://www.publish.csiro.au/an/an16577 ) 
weight_gain_error<-0.145 # Goopy et al. (2018, https://www.publish.csiro.au/an/an16577 ) 
digestible_energy_error<-0.1
crude_protein_error<-0.1

# farm-scale variables

Farms<-unique(Kenyaherds$QQNO)
  
for (f in Farms)
  {coln<-paste0("Farm_",f)
  # update feeding system 
  scenarios[which(scenarios$Variable == "feeding_system"),coln]<-
     median(Kenyaherds[which(Kenyaherds$QQNO==f),"system"])
  
  # update digestible energy (mean for all animals for now) (with error)
  scenarios[which(scenarios$Variable == "DE_perc"),coln][1:2]<-
     mean(Kenyaherds[which(Kenyaherds$QQNO==f),"dig_energy"])*
                           c(1-digestible_energy_error,1+digestible_energy_error)
  
  # update crude protein (mean for all animals for now) (with error)
  scenarios[which(scenarios$Variable == "perc_crude_protein_diet"),coln][1:2]<-
     mean(Kenyaherds[which(Kenyaherds$QQNO==f),"CP."])*
                           c(1-crude_protein_error,1+crude_protein_error)
  
  # update feed production CO2 (with error)
   scenarios[which(scenarios$Variable == "feed_prod_CO2"),coln][1:2]<-
     mean(Kenyaherds[which(Kenyaherds$QQNO==f),"feed_kgCO2"])*
                           c(1-feed_prod_error,1+feed_prod_error)
  # update feed transport CO2 (with error)
   scenarios[which(scenarios$Variable == "feed_trans_CO2"),coln][1:2]<-
     mean(Kenyaherds[which(Kenyaherds$QQNO==f),"feed.transport_kgCO2"])*
                           c(1-feed_trans_error,1+feed_trans_error)
  # update milk yield for each cow type (animal type 4-6) (with error)
   for(typ in 4:6)
     if(sum(Kenyaherds$QQNO==f&Kenyaherds$animal_type==typ)>0)
       {scenarios[which(scenarios$Variable==paste0("milk_yield_",typ)),coln][1:2]<-
         mean(Kenyaherds[which(Kenyaherds$QQNO==f&Kenyaherds$animal_type==typ),
                         "milk_yield"])*c(1-milk_yield_error,1+milk_yield_error)
       scenarios[which(scenarios$Variable==paste0("CP_rate_",typ)),coln][1:2]<-
         mean(Kenyaherds[which(Kenyaherds$QQNO==f&Kenyaherds$animal_type==typ),
                         "Cp"]/0.1) *c(1,1)
   
     }
  # update live weight and weight gain for each animal type (with error)
   for(typ in 1:12)
     if(sum(Kenyaherds$QQNO==f&Kenyaherds$animal_type==typ)>0)
       {scenarios[which(scenarios$Variable==paste0("live_weight_LW_",typ)),coln][1:2]<-
         mean(Kenyaherds[which(Kenyaherds$QQNO==f&Kenyaherds$animal_type==typ),
                         "live_weight_LW"])*c(1-live_weight_error,1+live_weight_error)
       scenarios[which(scenarios$Variable==paste0("weight_gain_WG_",typ)),coln][1:2]<-
         mean(Kenyaherds[which(Kenyaherds$QQNO==f&Kenyaherds$animal_type==typ),
                         "weight_gain_WG"])*c(1-weight_gain_error,1+weight_gain_error)
     }
   
   for(varia in variables_to_update_posnorm)
   if(sum(is.na(scenarios[which(scenarios$Variable==varia),coln][1:2]))==0)   
     if(!equals(scenarios[which(scenarios$Variable==varia),coln][1],
               scenarios[which(scenarios$Variable==varia),coln][2]))
        scenarios[which(scenarios$Variable==varia&
                   scenarios$param=="distribution"),coln]<-"posnorm" else
            scenarios[which(scenarios$Variable==varia&
                   scenarios$param=="distribution"),coln]<-"const"         
   
   }
  
write.csv(scenarios,"data/farm_scenarios.csv",row.names = FALSE)

Here’s the resulting file (only the first 3 farms are shown, and only the first 20 rows):

knitr::kable(read.table("data/farm_scenarios.csv", sep = ",", fileEncoding = "UTF-8-BOM",header=TRUE)[1:20,1:5])
Variable param Farm_1 Farm_2 Farm_3
N_animal_type_1 both 0 1 1
N_animal_type_2 both 0 0 0
N_animal_type_3 both 0 0 0
N_animal_type_4 both 0 1 0
N_animal_type_5 both 0 0 2
N_animal_type_6 both 2 0 0
N_animal_type_7 both 1 0 0
N_animal_type_8 both 0 0 0
N_animal_type_9 both 0 1 0
N_animal_type_10 both 0 0 0
N_animal_type_11 both 0 0 0
N_animal_type_12 both 0 0 0
feeding_system both 2 3 3
milk_yield_4 lower 2.826069862825
milk_yield_4 upper 4.969984931175
milk_yield_4 distribution posnorm
milk_yield_5 lower 3.8373469091625
milk_yield_5 upper 6.7484376678375
milk_yield_5 distribution posnorm
milk_yield_6 lower 5.471348045275

The first few rows in this table adjust the herd composition parameters, so that the scenario herd corresponds to the herd of the specific household. Other parameters including the milk yield per cow are also adjusted.

Now we can run the Monte Carlo simulation with these scenarios.

GHG_simulation_scenarios<-scenario_mc(base_estimate = estimate_read_csv("data/livestock_ghg_input_table.csv"),
                                 scenarios = read.csv("data/farm_scenarios.csv",
                                                      fileEncoding = "UTF-8-BOM"),
                                 model_function = ghg_emissions,
                                 numberOfModelRuns = 10,
                                 functionSyntax = "plainNames")

Given these results we can now plot distributions of the various outputs from the model, using the plot_distributions function.

# enteric_CH4=enteric_CH4,
enteric_CH4_dist_plot <- plot_distributions(mcSimulation_object = GHG_simulation_scenarios,
                   vars = "enteric_CH4",
                   method = "boxplot_density",
                   old_names = "enteric_CH4",
                   new_names = "CH4 Emissions from Enteric Fermentation per Farm") + 
  theme(axis.title.x=element_blank(),
        axis.title.y=element_blank())

# milk_yield=sum(total_milk),
milk_yield_dist_plot <- plot_distributions(mcSimulation_object = GHG_simulation_scenarios,
                   vars = "milk_yield",
                   method = "boxplot_density",
                   old_names = "milk_yield",
                   new_names = "Total Milk Yield") + 
  theme(axis.title.x=element_blank(),
        axis.title.y=element_blank())

# enteric_CH4_CO2eq=sum(enteric_CH4_CO2eq),
enteric_CH4_CO2eq_dist_plot <- plot_distributions(mcSimulation_object = GHG_simulation_scenarios,
                   vars = "enteric_CH4_CO2eq",
                   method = "boxplot_density",
                   old_names = "enteric_CH4_CO2eq",
                   new_names = "Carbon Dioxide Equivalent Enteric Methane Emissions from Livestock") + 
  theme(axis.title.x=element_blank())

# mm_CH4_CO2eq=sum(mm_CH4_CO2eq),
mm_CH4_CO2eq_dist_plot <- plot_distributions(mcSimulation_object = GHG_simulation_scenarios,
                   vars = "mm_CH4_CO2eq",
                   method = "boxplot_density",
                   old_names = "mm_CH4_CO2eq",
                   new_names = "Carbon Dioxide Equivalent Methane Emissions from Livestock Manure") + 
  theme(axis.title.x=element_blank(),
        axis.title.y=element_blank())

# mm_N2O_CO2eq=sum(mm_N2O_CO2eq),
mm_N2O_CO2eq_dist_plot <- plot_distributions(mcSimulation_object = GHG_simulation_scenarios,
                   vars = "mm_N2O_CO2eq",
                   method = "boxplot_density",
                   old_names = "mm_N2O_CO2eq",
                   new_names = "Carbon Dioxide Equivalent Nitrous Oxide Emissions from Livestock Manure") + 
  theme(axis.title.x=element_blank(),
        axis.title.y=element_blank())

# feed_CO2=feed_CO2,
feed_CO2_dist_plot <- plot_distributions(mcSimulation_object = GHG_simulation_scenarios,
                   vars = "feed_CO2",
                   method = "boxplot_density",
                   old_names = "feed_CO2",
                   new_names = "Carbon Dioxide Emissions from Livestock Feed Production and Transport") + 
  theme(axis.title.x=element_blank(),
        axis.title.y=element_blank())

# on_farm=sum(on_farm)
on_farm_dist_plot <- plot_distributions(mcSimulation_object = GHG_simulation_scenarios,
                   vars = "on_farm",
                   method = "boxplot_density",
                   old_names = "on_farm",
                   new_names = "Total Carbon Dioxide Equivalent Emissions from Livestock") + 
  theme(axis.title.y=element_blank())

library(patchwork)
## 
## Attaching package: 'patchwork'
## The following object is masked from 'package:raster':
## 
##     area
# remove the axis marks from the two lefthand-most plots and the x-axis title from the left and right plots

# plot with patchwork
layout <- "
           AAABBB
           CCCDDD
           EEEFFF
           GGGGGG
          "

enteric_CH4_dist_plot + milk_yield_dist_plot + 
  enteric_CH4_CO2eq_dist_plot + mm_CH4_CO2eq_dist_plot + 
  mm_N2O_CO2eq_dist_plot + feed_CO2_dist_plot + 
  on_farm_dist_plot + 
plot_layout(guides = "collect", design = layout) 

ggsave("Fig_1_GHG_simulation_outputs.png")
## Saving 10 x 10 in image

# THIS STILL NEEDS IMPROVEMENT - e.g. units are missing, and plot titles don’t fit (use concise symbols with subscripts etc.)

Possible errors resulting from input uncertainty

To illustrate the possible influence of input uncertainty on estimates of total farm emissions, we show the example of farm #5 from the household survey. This farm has 6 head of cattle, including 3 lactating cows. We simulate 1000 versions of this farm.

scenarios<-read.csv("data/farm_scenarios.csv", fileEncoding = "UTF-8-BOM")
farm_5<-scenarios[,c(1,2,7)]

farm5_emissions<-scenario_mc(base_estimate = estimate_read_csv("data/livestock_ghg_input_table.csv"),
                             scenarios = farm_5,
                             model_function = ghg_emissions,
                             numberOfModelRuns = 1000,
                             functionSyntax = "plainNames")

plot_distributions(mcSimulation_object = farm5_emissions,
                   vars = "on_farm",
                   method = "boxplot_density",
                   old_names = "on_farm",
                   new_names = "Total Carbon Dioxide Equivalent Emissions from Livestock") + 
  theme(axis.title.y=element_blank())

ggsave("Fig_2_uncertainty_illustration.png")
## Saving 7 x 5 in image

# NEED TO ADD THE USUAL PRECISE VALUE HERE FOR COMPARISON (also need to add units and possibly fix up other things)

This

farm5_emissions$y[,"CO2eq_per_milk"]<-farm5_emissions$y$pc_on_farm/farm5_emissions$y$pc_milk_yield

ggplot(data=farm5_emissions$y, aes(CO2eq_per_milk)) +
geom_histogram(aes(y=..density..), colour="black", fill="white")+
 geom_density(alpha=.2, fill="#FF6666") +
  ylab("Relative probability") +
    xlab(expression(Emission~density~of~milk~production~(kg~CO[2]-eq~kg^-1)))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggsave("Fig_3_emission_density_farm5.png")
## Saving 7 x 5 in image
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# STILL NEED TO add the usual precise value estimate

Population-wide emission density of milk production

Since the emission density of each household thus involves considerable uncertainty, such uncertainty should also be expected at the population level.

We can illustrate this uncertainty using a scatter plot. The following figure shows all combinations of total milk yield and total greenhouse gas emissions across all households.

As with the GHG Model, we apply the uncertainty functions on the milk_yield and on_farm variables. These come from the model results (4140) of the GHG_simulation_scenarios model results to the in_var (x) and an outcome variable out_var (y). In the model we defined on_farm as the sum of Carbon Dioxide Equivalent Enteric Methane Emissions from Livestock (enteric_CH4_CO2eq), Carbon Dioxide Equivalent Methane Emissions from Livestock Manure (mm_CH4_CO2eq), Carbon Dioxide Equivalent Nitrous Oxide Emissions from Livestock Manure (mm_N2O_CO2eq) and Carbon Dioxide Emissions from Livestock Feed Production and Transport (feed_CO2).

uncertainty::varscatter(in_var = GHG_simulation_scenarios$y$pc_milk_yield,
                        out_var = GHG_simulation_scenarios$y$pc_on_farm,
                        xlab = expression(Total~milk~yield~(kg~head^{-1}~year^{-1})),
                        ylab = expression(Total~emissions~(kg~CO[2]-eq~head^{-1}~year^{-1})))

## [1] "Scatter plot of influencing variable (x) and outcome variable (y) with associated and estimated densities."

# NICE PLOT, but would be nicer still if this was also a ggplot figure. Kinda inconsistent if it isn’t. Also can’t save it with ggsave, so will need another way.

Emissions intensity of milk production

An indicator that is often of interest is the ratio between these two variables: total emissions divided by the milk yield. We’ll calculate this next. We then have to remove all households that don’t produce any milk, because this ratio would be infinity for them. In fact, we’ll remove all households that reported very low milk yields of <500 kg per year per capita, since these are unlikely to represent regularly operating dairy operations. For instance, they may be mostly focused on beef cattle production, or they may only have a small number of dairy cows that happened to be unproductive at the time of the survey.

GHG_simulation_scenarios$y[,"CO2eq_per_milk"]<-GHG_simulation_scenarios$y$pc_on_farm/GHG_simulation_scenarios$y$pc_milk_yield

normal_milk_yield<-which(GHG_simulation_scenarios$y$pc_milk_yield>500)

dairy_households<-GHG_simulation_scenarios
dairy_households$x<-dairy_households$x[normal_milk_yield,]
dairy_households$y<-dairy_households$y[normal_milk_yield,]

Now we can plot the emissions intensity vs. the milk yield.

uncertainty::varscatter(in_var = dairy_households$y$pc_milk_yield,
                        out_var = dairy_households$y$CO2eq_per_milk,
                        xlab = expression(Total~milk~yield~(kg~head^{-1}~year^{-1})),
                        ylab = expression(Total~emissions~intensity~of~milk~production~(kg~CO[2]-eq~kg^{-1})))

## [1] "Scatter plot of influencing variable (x) and outcome variable (y) with associated and estimated densities."

This can also be illustrated by a density surface.

uncertainty::varkernel(in_var = dairy_households$y$pc_milk_yield,
                        out_var = dairy_households$y$CO2eq_per_milk,
                        xlab = expression(Total~milk~yield~(kg~head^{-1}~year^{-1})),
                        ylab = expression(Total~emissions~intensity~of~milk~production~(kg~CO[2]-eq~kg^{-1})))

## [1] "Density surface plot of an expected influential variable (x) and an outcome variable of interest (y)."

The milk yield per capita is sometimes proposed as an indicator of greenhouse gas emissions. We can now examine the merits of this indicator. Essentially, its usefulness hinges on its relationship with greenhouse gas emission intensity. We can investigate this by taking slices through the density surface.

Probability of GHG emissions for a given milk yield

The following plot shows the emissions intensity that corresponds to a milk yield of 2000 kg per capita per year, based on the population of households covered in the survey.

slice<-uncertainty::varslice_resample(in_var = dairy_households$y$pc_milk_yield,
                                 out_var = dairy_households$y$CO2eq_per_milk, 
                                 expectedin_var=2000,
                           n = 1000,
                           n_samples = 1000,
                           out_var_sampling = 1000) 

graphics::plot(x = slice$slice$Output_values, 
                 y = slice$slice$Relative_probability, 
                 type = "l", 
                 ylab = "Relative probability", 
                 xlab = "Emissions intensity", 
                 col = "seagreen", 
                 lwd = 2)

ggplot(data=data.frame(CO2eq_per_milk=slice$resampled), aes(CO2eq_per_milk)) +
geom_histogram(aes(y=..density..), colour="black", fill="white")+
 geom_density(alpha=.2, fill="#FF6666") +
  ylab("Relative probability") +
    xlab(expression(Emission~density~of~milk~production~at~2000~kg~milk~per~capita~(XX~CO[2]-eq~kg^-1)))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggsave("Fig_X_emission_density_all_farms.png")
## Saving 7 x 5 in image
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

We see considerable uncertainty about the emissions at this particular level of production, and the previous figure shows that emissions of all other levels is similarly variable. This raises considerable doubts about the usefulness of milk yield as an indicator of emissions intensity - at least when this indicator is taken at face value, without considering this uncertainty. Yet, some information is contained even in such an imperfect indicator, since, as also shown in the previous figure, there is clearly a correlation between milk yield and emissions intensity. The following proposes a mechanism for quantifying the level of certainty, with which a change in milk yield can reasonably be associated with a change in emissions intensity.

Probability of an increase in milk yield being associated with an increase in GHG emission intensity

To derive the level of certainty, with which a change in milk yield is associated with an increase or decrease in emissions intensity, we have to investigate the probability distributions (as shown in Fig. X) of emissions intensity at two levels of milk production. This can be achieved using the varslice_resample function for both production levels.

Let’s look at this for production of 2000 and 3000 kg milk per capita per year. We’ll take 1000 random draws of the emissions intensities that associated with both levels and compare values among these samples.

slice_2000<-uncertainty::varslice_resample(in_var = dairy_households$y$pc_milk_yield,
                                 out_var = dairy_households$y$CO2eq_per_milk, 
                                 expectedin_var=2000,
                           n = 1000,
                           n_samples = 1000,
                           out_var_sampling = 1000)
slice_3000<-uncertainty::varslice_resample(in_var = dairy_households$y$pc_milk_yield,
                                 out_var = dairy_households$y$CO2eq_per_milk, 
                                 expectedin_var=3000,
                           n = 1000,
                           n_samples = 1000,
                           out_var_sampling = 1000) 

sum(slice_3000$resampled < slice_2000$resampled) / length(slice_2000$resampled)
## [1] 0

The comparison in the last line of code above reveals that in 0 % of cases, emissions intensity at the higher milk production level is lower than at the lower level.

We can now generalize this result by doing similar calculations for multiple production levels and illustrating them for the whole population.

in_var = dairy_households$y$pc_milk_yield
out_var = dairy_households$y$CO2eq_per_milk


milk_yields_to_sample<-seq(100,10000,100)
milk_yields_to_sample<-milk_yields_to_sample[which(milk_yields_to_sample<max(in_var))]
milk_yields_to_sample<-milk_yields_to_sample[which(milk_yields_to_sample>min(in_var))]
slices<-list()

# slice through the distribution at every 100 L milk yield per cow

for(i in 1:length(milk_yields_to_sample))
{slices[[i]]<-
  varslice_resample(in_var,
                      out_var,
                      expectedin_var = milk_yields_to_sample[i],
                      n_samples = 100000,
  out_var_sampling = 1000)$resample
 }

decreased_emissions_probability <-
  data.frame(Base_yield = NA,new_yield=NA,probability=NA)

for (i in 1:(length(milk_yields_to_sample) - 1))
{
  for (j in (i + 1):length(milk_yields_to_sample))
    {decreased_emissions_probability[nrow(decreased_emissions_probability)+1,"Base_yield"]<-milk_yields_to_sample[i]
     decreased_emissions_probability$new_yield[nrow(decreased_emissions_probability)]<-milk_yields_to_sample[j]
     decreased_emissions_probability$probability[nrow(decreased_emissions_probability)] <-
      sum(slices[[i]] > slices[[j]]) / length(slices[[i]])
  
}
}

decreased_emissions_probability<-decreased_emissions_probability[which(!is.na(decreased_emissions_probability$Base_yield)),]

decreased_emissions_probability[,"yield_increase"]<-decreased_emissions_probability$new_yield-decreased_emissions_probability$Base_yield




ggplot(data = decreased_emissions_probability, aes(Base_yield, yield_increase, z =
                                                      probability)) + geom_contour_filled() +
  ggtitle("Probability that emission intensity reductions occurred") + 
  xlab("Baseline milk yield per head and year (kg)") +
  ylab("Milk yield increase per head and year (kg)") +
  stat_contour(breaks = c(0.9),col="black")
## Warning: stat_contour(): Zero contours were generated
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf

ggsave("Fig_X1_probability_of_intensity_reduction.png")
## Saving 7 x 5 in image
## Warning: stat_contour(): Zero contours were generated
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf

The 90% confidence contour is shown in the plot. The results indicate that whenever milk yields increase, the probability that emissions intensity decreased is usually greater than 50% (slightly lower values are due to randomness in the resampling procedure), but that large increases are needed to be sure of this.

Testing interventions - strategy

Farms may try to reduce their emissions intensity by making changes to their herd structure, feeding strategy or management practices. Programs aiming to effect emissions reductions in dairy production may want to reward such changes, and they may be looking to the milk production per capita as an indicator for the effectiveness of such measures. We now want explore whether this is a reasonable strategy.

Strategy

Develop scenarios that implement innovations

  • breed adjustment
  • herd structure

Quantify improvements on farms by comparing with the baseline. Plot improvements vs. the baseline on the above figure. Plot improvements vs. baseline in the emissions intensity figure.

For each slice, and each level of confidence, calculate the percentage of farms that made changes that would be correctly rewarded vs. the ones that did not change (the baseline population). Make a figure showing the Type I and Type II errors resulting from each confidence level at each baseline yield.

Testing interventions

Single-run function

To clearly characterize the impact of interventions, we should compare their impact with a counterfactual, i.e. with the same farms without the intervention. While this is usually impossible in practice, it is easy with a model. We just have to run the simulation again with identical inputs except for those that are affected by the intervention.

We already simulated emissions for a large population of farms, and we can now apply our interventions to this population by modifying certain parameters of the input tables. What we don’t want now, however, is to randomize all values again, since this would add additional noise to the results, which would mask the signal caused by the interventions.

We now first need a function that can apply our greenhouse gas emissions model to a pre-defined set of input values. We can find syntax that can achieve this in the mcSimulation function. Here it is (with slight modifications):

run_model<-function(x_data,model)
  {
  varnames<-colnames(x_data)[which(!colnames(x_data)=="Scenario")]
  model_function_ <- function(x) {
    
    e <- as.list(sapply(X = varnames, FUN = function(i) as.numeric(unlist(x[i]))))
    eval(expr = body(model), envir = e)
  }
  y <- do.call(what = rbind, args = lapply(X = apply(X = x_data, 
                                                     MARGIN = 1, FUN = model_function_), FUN = unlist))
  returnObject <- list(y = data.frame(y), x = data.frame(x_data))
  returnObject$call <- match.call()
  class(returnObject) <- cbind("mcSimulation", class(returnObject))
  return(returnObject)
}

The following code can precisely reproduce the results of the previous analysis:

res<-run_model(GHG_simulation_scenarios$x,ghg_emissions)

To simulate interventions, we can now feed this function with modified x_data data.frames.

Breed intervention

Introduction of higher-performance breeds may reduce greenhouse gas emissions. We simulate an intervention that replaces animals of low-performance breeds with animals of high-performance breeds.

To retain the farm-specific performance levels obtained from the farm survey, we simulate these impacts in a relative manner. For this, we first need to determine the performance and characteristics of each of the breeds. We’ll base this on the farm survey data (Kenyaherds dataset from above).

cows<-Kenyaherds[which(Kenyaherds$animal_type %in% 4:6),]

ggplot(cows, aes(factor(breed),milk_yield)) +
  geom_violin(draw_quantiles=c(0.5),fill="light green") +
  xlab("Breed") + ylab("Milk yield (kg per day)") + theme_bw()
## Warning: Groups with fewer than two data points have been dropped.

milk_summary<-cbind(aggregate(cows$milk_yield, by=list(cows$breed), FUN=median),table(cows$breed))
milk_summary<-milk_summary[,c(1,4,2)]
colnames(milk_summary)<-c("Breed","n","Median_milk_yield")
milk_summary
##   Breed   n Median_milk_yield
## 1     1 488          6.708090
## 2     2 134          4.999040
## 3     3  16          6.372821
## 4     4  20          5.538156
## 5     5  36          3.035786
## 6     7   1          7.217000
## 7     8   7          3.507483

This evaluation shows that the difference in median milk yields is relatively small, but that the milk yield distributions show a lot of variation that is apparently not explained by the breed. It is also apparent that the number of animals these estimates are based on varies considerably. Two breeds make up 88.6 % of the entire population, with others having very few animals, which means that the medians may be somewhat unreliable (especially for the single animal of breed 7).

Let’s define the intervention as follows:

All animals are replaced by high-performance animals of breed 1, which increases milk yields proportionally to the median milk yield differences. This means that milk yields are multiplied by the following factors:

milk_summary[,"Intervention_factor"]<-sapply(milk_summary$Median_milk_yield, function(x) min(x/milk_summary$Median_milk_yield[1],1))
                                         
milk_summary             
##   Breed   n Median_milk_yield Intervention_factor
## 1     1 488          6.708090           1.0000000
## 2     2 134          4.999040           0.7452256
## 3     3  16          6.372821           0.9500201
## 4     4  20          5.538156           0.8255936
## 5     5  36          3.035786           0.4525559
## 6     7   1          7.217000           1.0000000
## 7     8   7          3.507483           0.5228735

A change of the breed also causes changes in weight, mature weight and weight gain.

live_weights<-aggregate(Kenyaherds$live_weight_LW, by=list(Kenyaherds$animal_type,Kenyaherds$breed), FUN=median)
colnames(live_weights)<-c("animal_type","breed","median_live_weight")
live_weights[,"Intervention_factor"]<-NA
for(i in 1:nrow(live_weights))
  live_weights$Intervention_factor[i]<-live_weights$median_live_weight[i]/
     live_weights$median_live_weight[which(live_weights$animal_type==live_weights$animal_type[i] &                                                                                                            live_weights$breed==1)]

# mature weight should probably be changed, and the code below does this. However, the mature weight is currently the same for all the breeds, so this accomplished nothing.

mature_weights<-aggregate(Kenyaherds$mature_weight_MW, by=list(Kenyaherds$animal_type,Kenyaherds$breed), FUN=median)
colnames(mature_weights)<-c("animal_type","breed","median_mature_weight")
mature_weights[,"Intervention_factor"]<-NA
for(i in 1:nrow(mature_weights))
  mature_weights$Intervention_factor[i]<-mature_weights$median_mature_weight[i]/
     mature_weights$median_mature_weight[which(mature_weights$animal_type==mature_weights$animal_type[i] &                                                                                                            mature_weights$breed==1)]

# same for the weight gain - doesn't depend on the breed at the moment, so all the factors are 1

weight_gains<-aggregate(Kenyaherds$weight_gain_WG, by=list(Kenyaherds$animal_type,Kenyaherds$breed), FUN=median)
colnames(weight_gains)<-c("animal_type","breed","median_weight_gain")
weight_gains[,"Intervention_factor"]<-NA
for(i in 1:nrow(weight_gains))
  weight_gains$Intervention_factor[i]<-weight_gains$median_weight_gain[i]/
     weight_gains$median_weight_gain[which(weight_gains$animal_type==weight_gains$animal_type[i] &                                                                                                            weight_gains$breed==1)]
weight_gains$Intervention_factor[which(weight_gains$median_weight_gain==0)]<-1

Now we can modify the input table accordingly.

We go through all the farms and determine the ‘improvement’ coefficients (based on herd structure).

farms<-unique(Kenyaherds$QQNO)
breed_intervention<-data.frame(Farm=farms, milk_yield_4=1, milk_yield_5=1, milk_yield_6=1) 

Kenyaherds_interventions<-Kenyaherds[,c("QQNO",
                                        "system",
                                        "animal_type",
                                        "breed",
                                        "live_weight_LW",
                                        "mature_weight_MW",
                                        "weight_gain_WG",
                                        "milk_yield",
                                        "dig_energy",
                                        "CP.")]

Kenyaherds_interventions[,"milk_change"]<-0
cows<-which(Kenyaherds_interventions$animal_type %in% 4:6)
Kenyaherds_interventions$milk_change[cows]<- sapply(cows,function(x) 1/milk_summary$Intervention_factor[milk_summary$Breed==Kenyaherds_interventions$breed[x]]  )
  
milk_change<-aggregate(Kenyaherds_interventions$milk_change[cows], by=list(Kenyaherds_interventions$QQNO[cows]), FUN=mean)

breed_intervention$milk_yield_4[which(breed_intervention$Farm %in% milk_change[,1])]<-milk_change[,2]
breed_intervention$milk_yield_5<-breed_intervention$milk_yield_4
breed_intervention$milk_yield_6<-breed_intervention$milk_yield_4

for(i in 1:nrow(Kenyaherds_interventions))
  {Kenyaherds_interventions[i,"live_weight_change"]<-
    1/live_weights$Intervention_factor[which(live_weights$animal_type==Kenyaherds_interventions$animal_type[i]&
                                             live_weights$breed==Kenyaherds_interventions$breed[i])]
#   Kenyaherds_interventions[i,"mature_weight_change"]<-
#      1/mature_weights$Intervention_factor[which(mature_weights$animal_type==Kenyaherds_interventions$animal_type[i]&
#                                                 mature_weights$breed==Kenyaherds_interventions$breed[i])]
   Kenyaherds_interventions[i,"weight_gain_change"]<-
     1/weight_gains$Intervention_factor[which(weight_gains$animal_type==Kenyaherds_interventions$animal_type[i]&
                                              weight_gains$breed==Kenyaherds_interventions$breed[i])]
}
   
live_weight_agg<-aggregate(Kenyaherds_interventions$live_weight_change,
                           by=list(Kenyaherds_interventions$QQNO,Kenyaherds_interventions$animal_type), FUN=mean)
#mature_weight_agg<-aggregate(Kenyaherds_interventions$mature_weight_change,
#by=list(Kenyaherds_interventions$QQNO,Kenyaherds_interventions$animal_type), FUN=mean)
weight_gain_agg<-aggregate(Kenyaherds_interventions$weight_gain_change,
                           by=list(Kenyaherds_interventions$QQNO,Kenyaherds_interventions$animal_type), FUN=mean)

colnames(live_weight_agg)<-c("farm","animal_type","conv_factor")
#colnames(mature_weight_agg)<-c("farm","animal_type","conv_factor")
colnames(weight_gain_agg)<-c("farm","animal_type","conv_factor")

for(i in 1:nrow(live_weight_agg))
  breed_intervention[which(breed_intervention$Farm==live_weight_agg$farm[i]),
                     paste0("live_weight_LW_",live_weight_agg$animal_type[i])]<-
    live_weight_agg$conv_factor[i]

#for(i in 1:nrow(mature_weight_agg))
#  breed_intervention[which(breed_intervention$Farm==mature_weight_agg$farm[i]),
#                     paste0("mature_weight_LW_",mature_weight_agg$animal_type[i])]<-
#    mature_weight_agg$conv_factor[i]

for(i in 1:nrow(weight_gain_agg))
  breed_intervention[which(breed_intervention$Farm==weight_gain_agg$farm[i]),
                     paste0("weight_gain_WG_",weight_gain_agg$animal_type[i])]<-
    weight_gain_agg$conv_factor[i]

breed_intervention[is.na(breed_intervention)]<-1

Go through the x file, for each row take a random draw based on adoption rate (does farm adopt or not), then (possibly) apply factors. Actually, in this case we want to simulate full adoption, so we’ll set this to 1.

x_breed_intervention <- GHG_simulation_scenarios$x

adoption_rate <- 1

for (i in 1:nrow(x_breed_intervention))
{
  if (rbinom(1, 1, adoption_rate)) # if farmers adopt, make changes in x file
    x_breed_intervention[i, colnames(breed_intervention)[2:ncol(breed_intervention)]] <-
      x_breed_intervention[i, colnames(breed_intervention)[2:ncol(breed_intervention)]] *
      breed_intervention[which(paste0("Farm_", breed_intervention$Farm) ==
                                 x_breed_intervention$Scenario[i]),
                         2:ncol(breed_intervention)]
}
# run model for modified x file
breed_res <- run_model(x_breed_intervention, ghg_emissions)

Compare outcomes of breed intervention with baseline

breed_res$y[,"CO2eq_per_milk"]<-breed_res$y$pc_on_farm/breed_res$y$pc_milk_yield

breed_res$y<-breed_res$y[which(GHG_simulation_scenarios$y$pc_milk_yield>500),]
breed_res$x<-breed_res$x[which(GHG_simulation_scenarios$y$pc_milk_yield>500),]

impact_of_breed<-data.frame(baseline=dairy_households$y$pc_milk_yield,breed_milk=breed_res$y$pc_milk_yield,CO2_intensity_baseline= dairy_households$y$CO2eq_per_milk,CO2_intensity_breed=breed_res$y$CO2eq_per_milk)
impact_of_breed[,"milk_change"]<-round(impact_of_breed$breed_milk-impact_of_breed$baseline,2)
impact_of_breed[,"emissions_intensity_change"]<-round(impact_of_breed$CO2_intensity_breed-impact_of_breed$CO2_intensity_baseline ,3)

ggplot(data=impact_of_breed,aes(x=milk_change,y=emissions_intensity_change)) + geom_point() + theme_bw() + ylab("Change in emissions intensity (kg CO2-eq kg-1)") + xlab("Change in per capita milk yield (kg)")

ggsave("Fig_XX_breed_reductions.png")
## Saving 7 x 5 in image
ggplot(data = decreased_emissions_probability, aes(Base_yield, yield_increase)) + geom_contour_filled(aes(z=probability)) +
  stat_contour(breaks = c(0.9),aes(z=probability),col="black") +
  geom_point(data=impact_of_breed[which(impact_of_breed$milk_change>0),],aes(x=baseline,y=milk_change)) +
  ggtitle("Probability that emission intensity reductions occurred") + 
  xlab("Baseline milk yield per head and year (kg)") +
  ylab("Milk yield increase per head and year (kg)") 
## Warning: stat_contour(): Zero contours were generated
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf

ggsave("Fig_XX1_breed_probability_of_intensity_reduction.png")
## Saving 7 x 5 in image
## Warning: stat_contour(): Zero contours were generated
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf

This plot shows a lot of linear patterns that result from the application of fixed factors to convert the milk yield between cows of different breeds. Not so pretty, but maybe ok.

Herd structure

Retiring unproductive bulls

To simulate this scenario, we simply have to remove all animals of type 1 from the herds.

x_nobull_intervention <- GHG_simulation_scenarios$x

nobull_adoption_rate <- 1.0

for (i in 1:nrow(x_nobull_intervention))
{
  if (rbinom(1, 1, nobull_adoption_rate)) # if farmers adopt, make changes in x file
    x_nobull_intervention$N_animal_type_1[i]<-0
}
# run model for modified x file
nobull_res <- run_model(x_nobull_intervention, ghg_emissions)

Evaluate impact of “Retiring unproductive bulls” intervention

We’re assuming that all bulls >36 months old have been removed.

uncertainty::varscatter(in_var = nobull_res$y$pc_milk_yield,
                        out_var = nobull_res$y$pc_on_farm,
                        xlab = expression(Total~milk~yield~(kg~head^{-1}~year^{-1})),
                        ylab = expression(Total~emissions~(kg~CO[2]-eq~head^{-1}~year^{-1}))) 
## [1] "Scatter plot of influencing variable (x) and outcome variable (y) with associated and estimated densities."
uncertainty::varscatter(in_var = GHG_simulation_scenarios$y$pc_milk_yield,
                        out_var = GHG_simulation_scenarios$y$pc_on_farm,
                        xlab = expression(Total~milk~yield~(kg~head^{-1}~year^{-1})),
                        ylab = expression(Total~emissions~(kg~CO[2]-eq~head^{-1}~year^{-1})))
## [1] "Scatter plot of influencing variable (x) and outcome variable (y) with associated and estimated densities."
uncertainty::varkernel(in_var = nobull_res$y$pc_milk_yield,
                        out_var = nobull_res$y$pc_on_farm,
                        xlab = expression(Total~milk~yield~(kg~head^{-1}~year^{-1})),
                        ylab = expression(Total~emissions~(kg~CO[2]-eq~head^{-1}~year^{-1})))
## [1] "Density surface plot of an expected influential variable (x) and an outcome variable of interest (y)."
uncertainty::varkernel(in_var = GHG_simulation_scenarios$y$pc_milk_yield,
                        out_var = GHG_simulation_scenarios$y$pc_on_farm,
                        xlab = expression(Total~milk~yield~(kg~head^{-1}~year^{-1})),
                        ylab = expression(Total~emissions~(kg~CO[2]-eq~head^{-1}~year^{-1})))
## [1] "Density surface plot of an expected influential variable (x) and an outcome variable of interest (y)."
No bulls intervention (left) vs. baseline (right)No bulls intervention (left) vs. baseline (right)No bulls intervention (left) vs. baseline (right)No bulls intervention (left) vs. baseline (right)

No bulls intervention (left) vs. baseline (right)

Reducing unproductive replacement males

The idea is to raise the off-take rate of replacement males by 23% (not sure why this number). We’ll implement this by eliminating replacement males on all adopting farms, with a 23% adoption rate. The result should be quite similar.

x_lessReplace_intervention <- GHG_simulation_scenarios$x

lessReplace_adoption_rate <- 0.23

for (i in 1:nrow(x_lessReplace_intervention))
{
  if (rbinom(1, 1, lessReplace_adoption_rate)) # if farmers adopt, make changes in x file
    x_lessReplace_intervention$N_animal_type_3[i]<-0
}
# run model for modified x file
lessReplace_res <- run_model(x_lessReplace_intervention, ghg_emissions)

Compare outcomes of herd intervention with baseline

No bulls scenario

nobull_res$y[,"CO2eq_per_milk"]<-nobull_res$y$pc_on_farm/nobull_res$y$pc_milk_yield

nobull_res$y<-nobull_res$y[which(GHG_simulation_scenarios$y$pc_milk_yield>500),]
nobull_res$x<-nobull_res$x[which(GHG_simulation_scenarios$y$pc_milk_yield>500),]

impact_of_nobull<-data.frame(baseline=dairy_households$y$pc_milk_yield,breed_milk=nobull_res$y$pc_milk_yield,CO2_intensity_baseline= dairy_households$y$CO2eq_per_milk,CO2_intensity_breed=nobull_res$y$CO2eq_per_milk)
impact_of_nobull[,"milk_change"]<-round(impact_of_nobull$breed_milk-impact_of_nobull$baseline,2)
impact_of_nobull[,"emissions_intensity_change"]<-round(impact_of_nobull$CO2_intensity_breed-impact_of_nobull$CO2_intensity_baseline ,3)

ggplot(data=impact_of_nobull,aes(x=milk_change,y=emissions_intensity_change)) + geom_point() + theme_bw() + ylab("Change in emissions intensity (kg CO2-eq kg-1)") + xlab("Change in per capita milk yield (kg)")

ggplot(data = decreased_emissions_probability, aes(Base_yield, yield_increase)) + geom_contour_filled(aes(z=probability)) +
  stat_contour(breaks = c(0.9),aes(z=probability),col="black") +
  geom_point(data=impact_of_nobull[which(impact_of_nobull$milk_change>0),],aes(x=baseline,y=milk_change)) +
  ggtitle("Probability that emission intensity reductions occurred") + 
  xlab("Baseline milk yield per head and year (kg)") +
  ylab("Milk yield increase per head and year (kg)") 
## Warning: stat_contour(): Zero contours were generated
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf

ggsave("Fig_XX1_breed_probability_of_intensity_reduction.png")
## Saving 7 x 5 in image
## Warning: stat_contour(): Zero contours were generated
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
ggsave("Fig_XX_impact_of_nobull_reductions.png")
## Saving 7 x 5 in image
## Warning: stat_contour(): Zero contours were generated
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf

Fewer replacement males scenario

lessReplace_res$y[,"CO2eq_per_milk"]<-lessReplace_res$y$pc_on_farm/lessReplace_res$y$pc_milk_yield

lessReplace_res$y<-lessReplace_res$y[which(GHG_simulation_scenarios$y$pc_milk_yield>500),]
lessReplace_res$x<-lessReplace_res$x[which(GHG_simulation_scenarios$y$pc_milk_yield>500),]

impact_of_lessReplace<-data.frame(baseline=dairy_households$y$pc_milk_yield,breed_milk=lessReplace_res$y$pc_milk_yield,CO2_intensity_baseline= dairy_households$y$CO2eq_per_milk,CO2_intensity_breed=lessReplace_res$y$CO2eq_per_milk)
impact_of_lessReplace[,"milk_change"]<-round(impact_of_lessReplace$breed_milk-impact_of_lessReplace$baseline,2)
impact_of_lessReplace[,"emissions_intensity_change"]<-round(impact_of_lessReplace$CO2_intensity_breed-impact_of_lessReplace$CO2_intensity_baseline ,3)

ggplot(data=impact_of_lessReplace,aes(x=milk_change,y=emissions_intensity_change)) + geom_point() + theme_bw() + ylab("Change in emissions intensity (kg CO2-eq kg-1)") + xlab("Change in per capita milk yield (kg)")

ggplot(data = decreased_emissions_probability, aes(Base_yield, yield_increase)) + geom_contour_filled(aes(z=probability)) +
  stat_contour(breaks = c(0.9),aes(z=probability),col="black") +
  geom_point(data=impact_of_lessReplace[which(impact_of_lessReplace$milk_change>0),],aes(x=baseline,y=milk_change)) +
  ggtitle("Probability that emission intensity reductions occurred") + 
  xlab("Baseline milk yield per head and year (kg)") +
  ylab("Milk yield increase per head and year (kg)") 
## Warning: stat_contour(): Zero contours were generated
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf

ggsave("Fig_XX1_breed_probability_of_intensity_reduction.png")
## Saving 7 x 5 in image
## Warning: stat_contour(): Zero contours were generated
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
ggsave("Fig_XX_lessReplace_reductions.png")
## Saving 7 x 5 in image
## Warning: stat_contour(): Zero contours were generated
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf

References

Dong, Hongmin, Joe Mangino, Tim McAllister, Jerry Hatfield, Donald Johnson, Keith Lassey, Magda Aparecida de Lima, and Anna Romanovskaya. 2006. IPCC Guidelines for National Greenhouse Gas Inventories - Volume 4 Agriculture, Forestry and Other Land Use. https://www.ipcc-nggip.iges.or.jp/public/2006gl/pdf/4_Volume4/V4_10_Ch10_Livestock.pdf.
IPCC. 2021. Good Practice Guidance and Uncertainty Management in National Greenhouse Gas Inventories. Intergovernmental Panel on Climate Change (IPCC). https://www.ipcc-nggip.iges.or.jp/public/gp/english/.
Luedeling, Eike, Lutz Goehring, Katja Schiffers, Cory Whitney, and Eduardo Fernandez. 2022. decisionSupport: Quantitative Support of Decision Making Under Uncertainty. http://www.worldagroforestry.org/.
R Core Team. 2022. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
Wilkes, Andreas, Charles Odhong, Suzanne van Dijk, Simon Fraval, and Shimels Eshete Wassie. 2019. “Methods and Guidance to Support MRV of Livestock Emissions (Working Paper No. 285).” CGIAR Research Program on Climate Change, Agriculture; Food Security (CCAFS).
Wilkes, Andreas, Shimels Wassie, Charles Odhong’, Simon Fraval, and Suzanne van Dijk. 2020. “Variation in the Carbon Footprint of Milk Production on Smallholder Dairy Farms in Central Kenya.” Journal of Cleaner Production 265 (August): 121780. https://doi.org/10.1016/j.jclepro.2020.121780.